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We report on a novel extension of the recent phase-field crystal (PFC) method introduced in [Elder 
et ai, Phys. Rev. Lett., 88, 245701:f-4 (2002)], which incorporates elastic interactions as well as 
crystal plasticity and diffusive dynamics. In our model, elastic interactions are mediated through 
wave modes that propagate on time scales many orders of magnitude slower than atomic vibrations 
but still much faster than diffusive times scales. This allows us to preserve the quintessential 
advantage of the PFC model: the ability to simulate atomic-scale interactions and dynamics on time 
scales many orders of magnitude longer than characteristic vibrational time scales. We demonstrate 
the two different modes of propagation in our model and show that simulations of grain growth and 
elasto-plastic deformation are consistent with the microstructural properties of nanocrystals. 

PACS numbers: 46.15.-x, 61.82.Rx, 62.25.+g, 62.30.+d, 63.22.+m 



The deformation of a solid triggers processes which op- 
erate across several length and time scales. On long 
length and time scales its behavior can be described 
by a set of hydrodynamic equations Q, Q, which de- 
scribe, e.g., elastic deformation dynamics of the body. On 
atomic length (~ 10~ 10 m) and time (~ 10~ 13 s) scales, 
on the other hand, the dynamics can be captured by 
direct molecular dynamics (MD) simulations, which in- 
corporate local bonding information either through di- 
rect quantum-mechanical calculations or semi-empirical 
many-body potentials. While innovations in computing 
methods have greatly improved the efficiency of MD sim- 
ulations, standard atomistic computer simulations are 
still limited to fairly small system sizes (~ 10 9 atoms) 
and short times (~ 10~ 8 s). This limitation is most severe 
when developing simulation models to study the physics 
and mechanics of nanostructured materials, where the 
relevant length scales are atomic and time scales are 
mesoscopic. In this regime, the available numerical tools 
are rare. 

Progress towards alleviating this limitation has re- 
cently been made by the introduction of a new modeling 
paradigm known as the phase-field crystal (PFC) method 
[3j. This method introduces a local atomic mass density 
field p(r) in which atomic vibrations have been integrated 
out up to diffusive time scales. Dissipative dynamics are 
then constructed to govern the temporal evolution of p. 
Unfortunately, the original PFC model evolves mass den- 
sity only on diffusive time scales. In particular, it does 
not contain a mechanism for simulating elastic interac- 
tions, an important aspect for studying, for example, the 
deformation properties of nanocrystallinc solids. 

In this Letter, we introduce a modified phase-field crys- 
tal (MPFC) model that includes both diffusive dynamics 
and elastic interactions. This is achieved by exploiting 
the separation of time scales that exists between diffusive 
and elastic relaxation processes in solids. In particular, 
the MPFC model is constructed to transmit long wave- 



length density fluctuations with wave modes that propa- 
gate up to a time scale t w , after which the strain-relaxed 
density field continues to evolve according to diffusive dy- 
namics. The key feature of our approach is that the value 
of t w can be chosen to be much smaller than the charac- 
teristic time scale of diffusion and still much larger than 
1/u>d ~ lCP 13 s, where ojd denotes the Debye frequency. 

The phase field crystal methodology begins by intro- 
ducing an effective free-energy expanded to lowest order 
in the mass density p(r): 

F[p; T] = J (p/2[r + (q + V 2 )> + p 4 /4)d 2 x, (1) 

Here, r ~ (T — T m )/Lcp and T m , L and cp are, respec- 
tively, the melting temperature, the latent heat of fusion 
and specific heat at constant pressure of the pure mate- 
rial. Also, q D = 2tt /a, where a is the equilibrium lattice 
spacing. This free energy is identical to the one used in 
Ref. 0] , and gives rise to a phase diagram of coexisting 
liquid, solid, and striped phases, as shown in Fig. 1(a). 
In the solid phase, p is non-zero everywhere and spatially 
periodic on the atomic scale with hexagonal symmetry in 
two spatial dimensions. In the liquid phase, p takes on a 
constant value everywhere. In the PFC formalism lattice 
sites are always occupied and vacancy diffusion and topo- 
logical defects are represented via modulations of the lo- 
cal density amplitude and wavelength. Elastic constants 
can be determined by computing Cyjy = 6 2 F/SuijSuki, 
where «,j represents the strain of a particular deforma- 
tion state. 

In the original PFC model, the evolution of the mass 
density is given by 

dp/dt = a 2 V 2 {SF[p; T]/5p) . (2) 

where a is a constant. A severe limitation of the PFC 
model in Eq.J2J is that it only allows for diffusive den- 
sity relaxation. The model does not inherently contain a 
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FIG. 1: (a) Phase diagram indicating isothermal quench. 
The diagram is symmetric around po = 0. (b-top) Snap- 
shot in the evolution of polycrystalline solidification using the 
MPFC model. Grain boundaries are highlighted in white, (b- 
bottom) Zoom-in of 4 crystal grains and their orientations. 

suitable separation of times scales between phase trans- 
formation kinetics and the much more rapid ("instan- 
taneous" ) elastic relaxation. This precludes the study of 
phase transformation phenomena in the presence of com- 
plex mechanical deformations. It should be pointed out 
that while homogeneous deformations can be imposed 
through an affine transformation |5| , this method is inap- 
plicable in cases where non-homogeneous stress distribu- 
tions arise. As will be demonstrated below, these serious 
shortcomings of the original PFC model can be circum- 
vented in a way that allows us to preserve the quintessen- 
tial advantage of the PFC model, namely, the ability to 
simulate atomic-scale interactions and dynamics on time 
scales many orders of magnitude longer than molecular 
dynamics time scales. Most importantly, our modified 
model naturally incorporates "instantaneous" elastic in- 
teractions. 

We begin by introducing a modified PFC equation 
given by 

d 2 p dp 22 

_ +/ 3_ = a Vp + £ (3) 

where p = SF[p;T]/Sp, while (3 and a are phenomeno- 
logical constants. F is the free-energy introduced in 
Eq. ^ and £ is a Gaussian random variable with correla- 
tions satisfying (££') = k B Tfiq d -^ /X 2 V 2 S(x-x')S(t-t'). 
Henceforth, we will set £ = 0. Equation |21 is of the form 
of a damped wave equation, containing two propagating 
density modes at early time and one diffusive mode at 
late times. Specifically, the fast dynamics of the MPFC 



model are governed by the first term of Eq. @ , while the 
late time dynamics are governed by Eq.(|5J|. 

To elucidate the dynamics described by Eq. we 
performed a Floquet stability analysis. This was done by 
assuming a perturbation in the density of the form p p = 

Peq + Sp, where p eq = p + Y, n , m a n,me' G " ,m ' r , with p Q the 
average density, G n ,m = nx+ (n + 2m) j \fiy the triangu- 
lar reciprocal lattice vectors and a n , m their correspond- 
ing amplitudes. Here, dp = E n , ra ^mW e °' m ' ?+,<3 '' , | 
where Q is a perturbation wave vector and b n>m (t) 
the perturbation amplitude of mode (m,n). Substitut- 
ing p p into the model and expanding to linear order 
gives an equation for 6 n , m . The leading order mode 
satisfies 6o.o ~ e tu}t . The dispersion relation u>(Q) 
is given by lo(Q) = i(3/2 ± A(Q)/2, where A (Q) = 
V-/? 2 + 4a 2 Q 2 [3p 2 +r + (Q2 - q*)* + 9/8A 2 n J. Note 
that when 4a 2 Q 2 [3p 2 + r + (Q 2 - q 2 Q ) 2 + 9/8A 2 mm ] > 
j3 2 , the dispersion is approximately ui(Q) ~ ift/2 ± 
2aQ ^3p 2 + r+(Q 2 - q 2 ) 2 + 9/8A 2 mm = i/3/2 ± v eff Q. 
This dispersion describes a pair of waves that propagate 
undamped for time t w w 2(3~ x and distance L ~ v e fftw, 
after which they become effectively diffusive as in Ref. |3j . 
It is precisely these propagating modes which mediate 
elastic interactions in the model. Details of this calcula- 
tion as will be presented elsewhere 0. 

This analysis demonstrates that Eq.Q allows us to 
transmit disturbances across long distances using wave 
modes that propagate on a time scale t w . This allows 
all atomic positions to relax to a position close to their 
deformed equilibrium positions prior to any significant 
diffusion taking place. Most notably, the time scale t w 
can be set many orders of magnitude larger than charac- 
teristic vibrational time scales (~ 10 _13 s) but still signif- 
icantly faster than the scale on which diffusive processes 
occur. 

We now turn to the treatment of the fully non-linear 
evolution of Eq. J3J • The details of our numerical proce- 
dures are as follows. All simulations were conducted on 
a rectangular grid using periodic boundary conditions. 
Space was measured in units of the lattice constant a, 
while the grid size Ax, time step At and coefficients a, (3 
were chosen according to the particular application. Ex- 
ternal loads were applied to the boundary of our system 
by using a penalty function method. In this method, an 
additional term, of the form P = M(x, y, t)y/(p — Pbdy) 2 , 
is added to the free energy. This term couples the sample 
density p to an imposed periodic density field pbdy- The 
support of pbdy is the same as the support of the func- 
tion M(x, y, t) > 0, which defines the shape of the desired 
loading surface. The form of P thus couples some portion 
of the sample's density (e.g. near the sample boundaries) 
to the imposed boundary potential pbdy , which results in 
the sample's density field becoming slaved to the peaks of 
Pbdy as \M(x,y,t)\ — > oo. As the applied potential field 
is translated, the sample's density field, along the load- 
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ing surfaces, adiabatically follows the applied field. This 
specific form of P also assures that our penalty function 
does not alter the equilibrium phase diagram of the basic 
free energy functional F[p; T] defined above. 

We first simulated isothermal solidification using 
Eq.(J3J) by preparing the system in the liquid state and 
subsequently setting the temperature below the coexis- 
tence line in the phase diagram. In this case, the dy- 
namics gives rise to nucleation and growth of the solid 
phase in the presence of thermal fluctuations. To facil- 
itate nucleation, several nucleation sites were initiated 
in the metastable liquid phase in the form of random 
(Gaussian) fluctuations. During solidification, we found 
that the effect of the first term in Eq.© was negligible, 
and the growth rates and morphology were essentially 
indistinguishable from those using Eq. (J2J . Figure H il- 
lustrates growth and impingement of several nuclei in 
an undercooled melt. The simulation was started with 
the liquid of average density po — 0.285 and dimension- 
less temperature r = —0.25; other parameters were set 
to (Ax, At, a, 0) = (vr/8, 0.001, 15, 0.9). The measured 
grain boundary energies per unit length are consistent 
with the usual Read-Shockley form 0,0 

To demonstrate the presence of elastic relaxation 
modes in the MPFC model, we performed simulations 
of an effectively one dimensional single-crystal speci- 
men under uniaxial tension. The system was prepared 
in the coexistence region as given by the phase di- 
agram, and the solid sample was surrounded by liq- 
uid. Model parameters used were (r, if), Ax, At, a, f3) — 
(-0.4, 0.31, tt/8, 0.001, 15, 0.9). A tensile load applied to 
a semi-infinite continuum elastic bar can be theoretically 
modelled as array of coupled masses and springs along 
the x-axis, as illustrated in Fig. [21 When an atom at 
the boundary is displaced by an amount D± to the left, 
a tensile stress wave will propagate to the right. When 
atomic oscillations stop, a linear displacement distribu- 
tion, D(x) = Dix/L, will be established along the bar. 
Plots of displacement vs. position in the case of con- 
stant strain rate applied to the boundary atom are shown 
in Fig. |3 at three different times. Here, the displace- 
ments were extracted by a peak tracking method, where 
the locations of local maxima in p were tabulated af- 
ter each time step. The data clearly shows that the re- 
sponse of the system is consistent with elasticity theory. 

To make contact with the previous PFC formulation 
in Eq.(j2J |4J, we repeated the same simulations with a 
ten- fold increase in the damping parameter j3 = 9. The 
computed displacements, plotted in the inset of Fig. 
show that the response becomes viscoelastic as damp- 
ing is increased. Therefore, Eq.© alone does not ad- 
equately describe elastic responses in strained crystals 
at finite strain rates, while Eq.Q naturally incorporates 
such phenomena. We note that for simple modes of de- 
formation, Eq.© can still be used to model elastic relax- 
ation by uniformly adjusting all the atomic peaks of the 
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FIG. 2: Schematic illustration of the atomic locations before 
and after tension is applied in the MPFC model. 
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FIG. 3: The displacements along a one-dimensional sample in 
simple uniaxial tension at three different times (top). Linear 
profiles are consistent with linear elasticity theory. Inset: a 
ten-fold increase in /3 leads to visco-elastic behavior. 



density field p after a specified number of numerical time 
steps, which is equivalent to carrying out an affine trans- 
formation. This method is applicable in, e.g;, elucidating 
the glide dynamics of a single dislocation |5j]. However, 
this approach cannot be used to handle elastic relaxation 
in the case of complexgeometries, non-uniform stresses, 
and high strain rates U0, Hoi Hl| . 
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FIG. 4: A portion of the sample used to examine disloca- 
tion glide velocity. Parameters used: (r, ip, Ax, At, a, /3) = 
(-1, 0.49, tt/4, 0.001, 15,0.9). 
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FIG. 5: Two regimes of dislocation glide. For high strain 
rates we observed continuous glide, while at lower strain rate 
the dislocation set into a stick-slip motion. Inset: Dislocation 
glide velocity Vs. applied strain rate. 




22a 



Lattice constant - a 



0.02 0.04 0.06 0.08 0.1 0.12 E yy 30* ^. 

FIG. 6: Strain concentration in a double notched sample 
under a uniaxial tension. Left: A strain map of the center 
portion of the sample displayed at the bottom. Boundary 
atoms are highlighted in black. Right: Plot represents a 
strain profile from the center of the sample into the root of 
the notch. The solid line is a guide to the eye. 



We also examined the dynamics of individual disloca- 
tions. The set-up for these simulations are shown in Fig. 
0J Specifically, the top part of the crystal initially con- 
tains N atoms and the bottom part N+l. After the sam- 
ple equilibrated an edge dislocation formed and a con- 
stant shear strain rate was applied. The time-averaged 
dislocation glide velocity v was found to be a linear func- 
tion of the strain rate 7, consistent with classical disloca- 
tion theory. This theory predicts that v = A f/(pdb), where 
Pd is the dislocation density and b is the magnitude of the 
Burger's vector [l^ . 

In order to elucidate the local dynamics of individual 
dislocations, we computed the average strain in the crys- 
tal as a function of time for different strain rates. These 
results, shown in Fig. |5J revealed two regimes of dislo- 
cation glide. The first was characterized by continuous 



glide (observed at large 7) and the second by a stick- 
slip gliding of the dislocation at low 7. In both cases 
the applied plastic strain was relieved by the motion of 
the dislocation, and the time-averaged strain remained 
constant. 

To further illustrate the properties of our MPFC 
model, the effect of uniaxial tension in a notched sam- 
ple was examined. Figure El shows that strain (stress) 
in a notched sample appropriately concentrates near 
the notches, as expected from linear elasticity theory. 
In particular, treating the case of a double notched 
plate the stress concentration for this geometry is K t — 
a yy aX / a yy ~ ^ 0> which is in excellent agreement with 
our simulation result 1.81 . It is noteworthy that a sim- 
ulation with the PFC model (Eq. GJ for the same sys- 
tem and using an affine transformation to approximate 
the strains in the sample, failed to produce the expected 
strain concentration. 

In conclusion, we have introduced a novel phase- 
field crystal model (MPFC), which extends the previous 
phase-field crystal formalism by generating dynamics on 
two time scales. Atomic positions are relaxed rapidly at 
early times in a manner consistent with elasticity theory, 
while late time dynamics are governed by diffusive dy- 
namics characteristic of phase transformation kinetics, 
vacancy diffusion, grain boundary kinetics and disloca- 
tion climb. It is expected that the MPFC model will 
help open a new window into the study of phase trans- 
formation kinetics and microstructure heterogeneity in 
high strain rate loading of nanocrystalline solids. 
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